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BACKGROUND OF THE INVENTION 

The calculation of sensitivity partial derivatives is a frequently occurring task during 
the engineering design and control process. Sensitivity models are derived from an 

1 5 underlying math model of the physical system. Traditionally there have been two approaches 
for developing the sensitivity partial derivatives: (i) analytical models, and (ii) numerical 
models using finite difference-like techniques. Analytical models are preferred when 
available, though the additional modeling, software development, and checkout are often 
time-consuming and expensive. Numerical methods are conceptually straightforward, 

20 however, many methods must sample a function two or more times to estimate the 
sensitivity. 

BRIEF DESCRIPTION OF THE DRAWING 

FIG. 1 illustrates a plot of relative growth rates of polynomials according to one 
25 example embodiment of the invention. 

FIG. 2 illustrates an example FOTRAN subroutine according to one example 
embodiment of the invention. 

FIG. 3 illustrates an example computer platform with an Object-Oriented Cartesian 
Embedding Algorithm (OCEA) according to one example embodiment of the invention. 

30 
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DETAILED DESCRIPTION OF THE INVENTION 



L INTRODUCTION 

The present invention presents a third alternative approach for generating the 
5 sensitivity data: An Object-Oriented Cartesian Embedding Algorithm (OCEA). The 

computational advantage of the OCEA approach is that a single function evaluation yields 
the function value and the associated Jacobian and Hessian partial derivatives. The user is 
completely hidden from the details of how the partial derivatives are generated. The 
algorithm has been implemented in Fortran 90 and Macsyma 2.4. Module functions (Fortran 

10 90/95) are used to encapsulate new data types, and extended math and library functions for 
handling vector, tensor, and embedded variables. According to certain example 
embodiments, the present invention supports math models for scalar, vector, linear matrix 
equations, and matrix inversion. 

OCEA is a chain rule-based evaluation technique for analytically evaluating partial 

1 5 derivatives with respect to input variables of functions defined by a high-level language 
[Turner, 2002]. All of the problem functions are assumed to have continuous partial 
derivatives through second-order. Operationally, OCEA replaces conventional computer 
algebra operations and functions with generalized OCEA versions of these capabilities. The 
use of operator-overloading techniques permits a familiar conceptual framework to surround 

20 the use of OCEA software. The operator-overloading techniques redefine the intrinsic 

mathematical binary operators { +, -, *, **, /} and unary functions { cos(x), sin(x), tan(x), 
log(x),. . ..}. New derived data types, interface operators, and generalized mathematical 
operators are developed for computing the first and second order partials. 

Automatic differentiation has existed as a research topic since the 1980's. The 

25 common idea has been to pre-process a programmed function, identify the unary and binary 
functional operations, and build up forward and backward computational strategies. A 
successful implementation of this approach is ADIFOR (Automatic Differentiation of 
FORtran), which transforms the model source code based on a specification of dependent and 
independent variables, and produces source code that calculates the derivatives of the 

30 dependent variables with respect to the independent variables [Griewank, 89; Bischof, Carle, 
Corliss, Griewank, and Hovland, 92; Bischof, Carle, Khademi, Mauer, and Hovland, 95; 
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Eberhard and Bischof, 96; Hovland and Heath 97]. Other related codes include AD01 [Pryce 
and Reid, 98] which handles higher order derivatives using the chain rule and computational 
graphs. OCEA, unlike these previously developed tools, has the capability to build models 
on the fly, exploit high order models, and provide simple extensions for handling matrix 

5 calculations. No distinction is made between the independent and dependent variables. 
OCEA operates at the elementary scalar level, regardless of the complexity of the 
mathematical object. 

Mathematically, OCEA replaces each scalar operation with a higher dimension 
object. For example, assuming the g is a lxl scalar in traditional computer algebra, and 

10 introducing the OCEA-based transformation process, one obtains: 

d . a 



g:=[g, Vg, VVg]; V = * — +.•• + ». 



dx x dx m ^ 



lxl mx\ 



where Xi (i = 1 , . . . ,m) denotes the vector of independent variables, h. = (S n , • * • , S im ) denotes a 

unit vector in the i-th coordinate direction, and the transformed version of g is now of 
dimension lx(l+m+m 2 ). Generalizations for higher dimensional versions of OCEA are 
15 obvious. Because of the rapid increase in the number of partial derivatives as the order of the 
OCEA method increases, it becomes critically important to exploit symmetry. As a point of 
comparison, if a q-th order version of OCEA is considered, the dimension of each 

transformed OCEA scalar becomes ^ m n = ( 1 - m q+l )/(l - m) when symmetry is not 

77=0 

exploited. If symmetry is exploited the dimension of each transformed OCEA scalar 



20 becomes ^ 



n=0\ W-l j 



\ 171 J 



where (*) denotes the binomial coefficient. A plot of the 



relative growth rates is presented in Figure 1 . 

A benefit of the OCEA operator overloading strategy is that standard programming 
language constructs are used in building mathematical models. The complier recognizes that 
new data types and automatically replaces the conventional intrinsic operators and functions 
25 with generalized OCEA intrinsic operators and functions. For example, imagine that one 
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wants to evaluate f = xyl4z using standard Fortran and OCEA-enhanced Fortran, as shown 
in Table 1 : 



Math Model 


/ = xy/yfz 


Fortran 

f := (scalar) 


f = x*yl sqrt(z) 


OCEA-Enhanced Fortran 

F := (scalar, vector, tensor) 


/ = • 


x*yl sqrt(z) 

A 

> a/ a/1 

dx dy dz 

~d 2 f d 2 / a 2 / 1 

dxdx dxdy dzdz 




Table 1: Comparison oi 


F Fortran and OCEA Fortran Models 
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where three independent OCEA variables are defined (namely, x,y,z) and the Fortran and 
OCEA-Enhanced Fortran models are seen to be identical. The OCEA calculations, however, 
are very different from the Fortran calculations. The analyst never formulates or codes the 
partial derivatives. OCEA represents a linearized Clifford algebra. 

10 This patent includes seven sections. Section 2 presents the mathematical foundation 

for OCEA. Section 3 presents the software operator overloading, data typing, interface 
operators, and generalized transformational operators for multiplication, division, and 
composite function evaluations. A differential equation sensitivity approach is presented in 
Section 4. A simple example is worked out in detail to highlight the operational steps 

1 5 required for using OCEA methodology in Section 5. Several non-trivial applications are 

presented in Section 6 for handling nonlinear vector functions, equation of motion generation 
using Lagrange's equations, and differential equations. Applications are presented in Section 
7. 
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II. MATHEMATICAL FORMULATION 

The mathematical models for the OCEA are presented in this section. An essential 
step in the algorithm is that the independent coordinates are transformed using Eq. (1), 
leading to: 











"0 • 


• 0" 
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0 • 


• 0 





5 lxl IjcI nxl nxn 

v v * 

ljc( l + /i + w 2 ) 

where S tj denotes the standard kronecker delta function and i = 1 . .,n. The independent 
variables initialize the calculations for the embedded variables. 



10 IL1 Intrinsic Operators and Functions 

The generalizations for the intrinsic mathematical operators and functions are 
presented for addition, subtraction, multiplication, division, and composite functions. 
Addition and subtraction are straightforward. The most challenging generalizations are 
required for the product, division, and composite function calculations. All of the following 

15 transform equations assume that Jacobian and Hessian data for g and h has been built up 
from previous computational operations. One should observe that the tensor parts of the 
operators above frequently require vector outer products to complete the mathematical 
models. Many opportunities exist for exploiting the sparse structure of the resulting vector 
and tensor operations. The embedded coordinate mathematical operators follow as: 



20 



25 



Addition: 



Subtraction: 



g + h = ^g, Vg, V Vg] + [h 9 VA, V Vh] 

= [g + K Vg + v/*, v Vg + v v/*] 

g - h - [g, Vg, VVg] - [K VA, V Vh] 
= [g-A,Vg-VA,VVg-VVA] 
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Product Rule: 



g*h :=[g, Vg, Wg]*[/*, Vh, VVh] 

= [g*h, g*Vh + h*Vg, g*VVh + Vg{vhf + Vh[v g J + h*VVgj 



Division Rule: 



glh := [g, Vg, Wg]/[/z, Vh, VVh] 



glh, 
r hVg-gVh \ 
*» J' 



/ - - 



/*VVg-gVV/?-Vg(v/j) - V/»(Vg) 



+ 



V 



10 



Composite Function Rule: 



/([«.%^«])j= 



/(«). f + 



15 where the three scalar values for f(g), df /dg , and d 2 f /dg 2 are computed using standard 
intrinsic mathematical functions. 



IL2 Example Embedded Function Calculations 

This section considers the calculation of embedded versions of <Jg and tan~ ! (g). 
20 The square root calculation is considered first. 
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II. 2. 1 yfg Calculation: Referring to the composite function rule above one has / (g) = , 

which leads to — = —== , and — — = — , so that the embedded composite square root 

dg yjg dg 2*g 

function is given by: 

Jg'-=yj[g,Vg, VVg] = 



r vg -vg(vg) VVg 



The tan (g) calculation is considered next. 

77.22 tan _1 (g) Calculation: Referring to the composite function rule above one has 

1 f 2 

10 f(g)~ tan "' (g) > which leads to — = , and — ^ = > so tnat tne embedded 

5g l+g 2 dg 2 (l + g 2 ) 2 

composite tan"'(g) function is given by: 



tan-'(g) := tan' 1 ([g, Vg, VVg]) = 



Similar generalizations have been developed for the operators "**", sin, cos, tan, asin, acos, 
15 exp, log, sinh, cosh, and tanh. Others can be easily developed as required. 

III. SOFTWARE ARCHITECTURE 

According to one example embodiment, the OCEA mathematical operators and 
functions are translated into an object-oriented language where operator-overloading 

20 techniques are employed to re-define the algebraic operations that control the computational 
procedures. According to one embodiment, A Fortran 90 version of OCEA may be 
provided. An OCEA program is shown executing on a suitable computing platform, such as 
a personal computer or workstation, or other computing system, in Figure 3, wherein the 
software is stored in a machine readable form in memory, on a magnetic or optical disk, or in 

25 the form of a digital signal in electrical form. Defined data types, interface operator 

definitions, and module functions are used for encapsulating the data types and providing a 
language extension that supports the automatic generation of partial derivatives. As shown in 
Section II, according to one example embodiment, three derived data types are required: 
vector, tensor, and embedded. The embedded capabilities provide the OCEA capabilities for 
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libraries of mathematical operators and functions. The introduction of derived data types 
provides a software advantage; namely, that the compiler can detect the data types involved 
in a calculation and invoke the correct subroutine or function without user intervention. This 
is achieved in the software architecture design by providing typed data subroutines and 
5 functions for covering all possible mathematical operations, (i.e., real and embedded, 
embedded and embedded, and so on). For example, Table 2 provides a partial list of the 
capabilities required for generalizing the assignment (=) and math operators. An additional 
advantage is that the compiler can use array-processing techniques for handling large 
problems. 
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Data Types (A&B) 


Operations 


Math Model 


Embedded, Embedded 


5 ~ 5 


A+B, A-B 


Embedded, Vector 




A=B 


Vector, Embedded 




A=B 


Vector, Scalar 


Itjjc" » « 

5 


A=b*A, A f =b, i=l,...,n 


Tensor, Tensor 




A=B 


Embedded, Embedded 


"sin" 


A=sin(B) 



Table 2: Operator Overloading for Data Type Operations 



Two detailed examples are provided for clarifying the procedure. The first example presents 
vector handling for adding two vectors and the second supports multiplying a tensor by a 

1 5 scalar. PV denotes the defined vector data type, where VP ART is the vector structure 

constructor variable. PT denotes the defined tensor data type, where TP ART is the tensor 
structure constructor variable. The variable ":" denotes an array assignment where all 
components are assigned. Mixed data types are not allowed. Explicit data typing is required 
for all variables. As a preprocessing step, the compiler finds the correct function by (i) 

20 identifying the operator, and (ii) checking for module procedures with the same the input 
data types. 

IIL1 Adding Vector Data Types 
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If the compiler detects that two vector data types are added and that the result is a 
vector data type, then an automatic link is generated those calls the function routine 
FUNCTION ADD _PV in Module PVHANDLING. Symbolically this equation looks like 

Detected by Interface Operator 

v v ' 

ft 

PVTYPE = PVTYPE + PV_TYPE 

u u 

Detected by Function Input 
Data Types 

> v ' 

5 

When vector data types are added, an interface operator for addition is defined and a name is 
assigned for the module procedure that takes two vector data types as input, as presented in 
the software fragment: 

10 MODULE PV_HANDLING ! Vector Operations 



INTEFACE OPERATOR (+) 

MODULE PROCEDURE ADD PV 
15 END INTERFACE 



FUNCTION ADD_PV(A,B) 

TYPE(PV) ADD_PV 
20 TYPE(PV), INTENT(IN)::A,B 

ADD_PV%VPART(:) = A%VPART(:) + B%VPART(:) 
END FUNCTION ADD PV 



25 END MODULE PV HANDLING 

IIL2 Multiplying a Tensor Data Type by a Real Scalar 

If the compiler detects that a real scalar multiplies a tensor data type and that the 
30 result is a tensor data type, an automatic link is generated that calls the function routine 
FUNCTION_R_PT in Module PT HANDLING. Symbolically this equation looks like 
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Detected by Interface Operator 

v v ' 

ft 

PTTYPE = SCALAR*PT_TYPE 

Detected by Function Input 
Data Types 



When a real variable and a tensor are multiplied, an interface operator is defined for 
multiplication and a name is assigned for the module procedure available for multiplying a 
5 tensor by a scalar, as presented in the software fragment: 



MODULE PTHANDLING ITensor Operations 



1 0 INTERFACE OPERATOR (*) 

MODULE PROCEDURE MULTRPT 
END INTERFACE 



15 FUNCTION MULT_R_PT(R, A) 

TYPE(PT) MULTRPT 
REAL(KIND=8), INTENT(IN)::R 
TYPE(PT), INTENT(IN)::A 
INTEGER: :I 

20 DOI=l,N 

MULT_R_PT%TPART(:,I) = R*A%TPART(:,I) 
END DO 
END FUNCTION MULT R PT 



END MODULE PT HANDLING 

The object oriented data typing hides the building of links to the specific subroutines and 
functions required for completing calculations. 
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IV. DIFFERENTIAL EQUATION SENSITIVITY 

Linked mechanical systems lead to first-order differential equations of the form: 

x = JI/(x)- ! /(*,*) (2) 
Partial derivatives of Eq. (2) are required for computing state transition matrices and high- 
dimensional sensitivity operators. Even at second-order the partial derivatives can become 
quite involved. Using OCEA versions of x, f, and M, one obtains 

x:=[M~7, v(Ar'/) s VV(M"7)] 

The value of the implicit solution is that the inverse of M(x) is never formed explicitly or 
manipulated to generate the partial derivatives. The complexity of the implicit solution is 
compared with a conventional sensitivity approach. 

IV A OCEA Gaussian Elimination Approach: 

The calculation for Eq. (2) is simplified by redefining the equation as 

M(x)x = f(x 9 t) (3) 

i 

where OCEA versions of the variables are defined as 

x:=[x,[x X) - x Xn ],[x XitXi x M x x ^]] 

f,x^x 2 "' ]] 

and 

M:=[M,[M, Xi - M^]\M^ M M - M^]] 

Equation (3) is solved by using an OCEA version of Gaussian elimination. The Gaussian 
elimination algorithm is identical to standard algorithms, except that a single Use 
EB_Handling Module command is introduced and variables are typed as needed. The 
generalized intrinsic operators handle all of the details. The resulting calculation builds the 
partial derivatives of the state differential equation, without having to form M" 1 or any matrix 
products. A numerical example is given in Section VI. 
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IV.2 Conventional Second Order Calculation: 



Analytically computing the second-order partial derivatives of Eq. (3) leads to 



5 



(4) 



where 




-i 



and 




-i 



10 By counting the operations for the second-order inverse matrix partials, one finds that 
(5n 2 + 1 In) 1 2 nxn matrix products must be formed (fully accounting for symmetry). 
Consequently, for a medium-sized system where n = 100, one must evaluate 25,550 100x100 
matrix products. Comparing Eqs. (2) and (4), one observes that the OCEA solution generates 
the right-hand-side of Eq. (4) without ever explicitly forming any of the complicated matrix 

15 products presented above. 

IK 3 State and Parameter Transition Matrix Calculations 

The state and parameter transition matrix calculations are presented in this section. 
Both analytical and OCEA models are presented for generating the partial derivatives. 
20 Because the transition matrix variables are implicitly defined, OCEA is used to generate the 
information required for building the math models for integrating the partials. For simplicity, 
only a first order version of OCEA is assumed for the derivations presented in Sections 
IV.3.1 through IV.3.3. 

25 IV J J Analytic State Transition Matrix: The state transition matrix is obtained by integrating 
a first order model of the form 



x = F(x,t), x Q =x[ 



(5) 
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where x is the nxl state vector and the sensitivity of the terminal state with respect to the 
initial state is required. Analytically the required partial derivatives are obtained by 
integrating Eq. (5) as: 



i 

X (t) = x(t 0 ) + \F(x(T),T)dT 



and computing the partial derivative with respect to x(to), leading to: 



dx(t) 
dx(t 0 ) 



dF dx(r) 
dx(r) dx(t 0 ) 



dr 



where the initial condition is defined by an nxn identity matrix. The governing differential 
equation for the state transition matrix is obtained by differentiating the equation above with 
respect to t, yielding the nxn matrix differential equation: 



d_ 
dt 



dF dx(t) . dx(t) 



dx(t) 

[dx(t 0 )) dx(t)dx(t o y dx(t 0 ) 



= 1 



(6) 



IV. 3. 2 Analytic Parameter State Transition Matrix: The parameter transition matrix is 
obtained by integrating a first order model of the form 

15 

x = F (x, t; p) , x a = x\ /=i , p is a (/nxl) vector of parameters (7) 

where x is the nxl state vector and the sensitivity of the terminal state with respect to the 
problem parameters is required. Analytically the required partial derivatives are obtained by 
integrating Eq. (7) as: 

20 x(t) = x(t 0 )+\F(x(T),r;p)dt 

'o 

and computing the partial derivative with respect to p, leading to: 

dp r 0 U*( r ) & d PJ 

where the initial condition is defined by an nxm zero matrix. The governing differential 
equation for the parameter transition matrix is obtained by differentiating the equation above 
25 with respect to t, yielding the nxm matrix differential equation: 
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d_ 

dt 



( dx^ 
\ d Pj 



dF dx(r) dF dx 



(8) 



dx(r) dp dp' dp 1 ' 0 



IV. 3. 3 OCEA State and Parameter Transition Matrix Calculations: The state transition 
5 matrix differential equation of Eq. (6) is obtained from an OCEA version of Eq. (5), as 
follows: 

x = F 

x := I^jc Vx J = [x, x 2 ] 

F:=j> VF] = [/J F 2 ] 

where the subscripts denote the different parts of the OCEA variable. By defining the partial 
derivative variable dx/dx 0 , the required differential equations are given by 

Xj = F { , x, | /q = x 0 

10 d dx „ 8x dx\ T ( 9 ) 



3x 0 5x 0 3x ( 



J r 0 «.tn 



0 



The parameter transition matrix differential equation of Eq. (8) is obtained from an OCEA 
version of Eq.(7), as follows: 
x = F 

x:=[x Vx ] = [* V x x V p x ]=[xj x 2 x 3 ] 
F:=[ F VF] = [F V x F V„F] = tf F * F A 
1 5 where OCEA variables now include both x and p— which allows the gradient operator to be 
split into state and parameter parts, where the subscripts denote the different parts of the 
OCEA variable. By defining the partial derivative variable dx/dp , the required differential 
equations are given by 



ddx_ dx dxs 00) 

~ F 2— + F 3 >—L -0 - 



dt dp dp dp 

20 For both Eqs. (9) and (10) OCEA builds the complex and tedious part of the calculation 
automatically. 
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V. SIMPLE EXAMPLE PROBLEM: 

To highlight the required embedded calculations, the following simplified 3D 
example application is presented. The goal is to compute an embedded version of the 

composite function y4xz , where the embedded versions of x, y, and z follow as 





f 




"0 


0 


0" 










x, [1,0,0], 


0 


0 


0 












0 


0 


0 






( x ) 








"0 


0 


0" 






y 






y, [0,1,0], 


0 


0 


0 














0 


0 


0 












"0 


0 


0" 










z, [0,0,1], 


0 


0 


0 












0 


0 


0 




) 



The following three steps are required: 

Step 1: Form the product x*z (using the product operator): 







"f 




"o 


0 


0" 








0 




0 


0 


0" 




x*z := 


X, 


0 




0 


0 


0 




* 




0 


> 


0 


0 


0 








_0_ 




0 


0 


0 








1 




0 


0 


0 







"0" 




"f 


x*z, X* 


0 


+ z* 


0 




1 




0 





"0 0 0" 




"1" 




"0" 


T 


"0" 




"f 


T 


"0 0 0" 


X* 


0 0 0 


+ 


0 


* 


0 


+ 


0 


* 


0 


+ z* 


0 0 0 




0 0 0 




0 




1 




1 




_0_ 




0 0 0_ 







z 




"o o r 






x*z, 


0 




0 0 0 








X 




10 0 
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Step 2: Compute the vx*z (using the composite function operator): 



Vx*~z := 



J. 



x*z, 



2*J. 



x*z 



-1 



4*(x*z 



.3/2 



2*Jx> 



0 0 1 

0 0 0 

1 0 0 




-z 



4*(x*z) 
0 



3/2 



0 



1 



x*z 



1 



x*z 



2*4x~*z~ 4*(x*zf 2 



2*4x*z 4*(x*z) 
0 0 
-x 2 

0 



3/2 



4*(x*z) 3/2 



Step 3: Compute the product y *Vx*z (using the product operator): 



y*\Jx*z 



: Vx*z, 



-y*z 




4*x*Vx*z 2*V**z 

2 , 0 



X*Z 



*Vx* 



4*x*z 
x 



2*>/ 
y*\fx*z 

4*x*z 2*Vx*z 4*z*V! 



Jt*Z 



x^z 



The calculations can become extremely complex; nevertheless, the user is completely hidden 
from the details of the calculations and the calculations are accurate to the working precision 
of the machine. 
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VL COMPLEX APPLICATIONS 



Three classes of applications have been tested using a Fortran 90 version of OCEA. 
The three application classes consist of: (1) Non-linear n-dimensional vector Jacobian and 
Hessian calculations. (2) Generation of Lagrange's equations of motion, and (3) Generation 
of state transition matrix partials. These examples are presented in Section VI. 1 through 
VI.3. 
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15 
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VI. 1 Nonlinear Function Evaluation 

A highly nonlinear 7x1 vector algebraic function, , of five independent 

variables, £ , is considered where F(%) is given by: 

e u +x 2 *cos(z) 

x*w*(jy*z) 2 
x*y*z 2 



F(4) = 



z 3 /u 



'yfw 



x*u*sin l (y/(u*w)} 
z 1/3 *ln(V^) 

and £ = (x, jy,z,w, w) . The computational procedure is outlined in Figure 2, which consists 

of a Fortran 90 subroutine for executing the OCEA calculations for the function defined 
above. The user inputs a real-valued vector of independent variable (i.e., INDVAR). These 
variables are internally transformed to OCEA variables for the calculation (i.e., EBJVAR). 
NV denotes the number of independent variables and NF denotes the number of function 
dimensions to be evaluated, where both are defined in Module Problem_Data. The Module 
EBJHandling provides all of the OCEA operator overloading and extended function 
capabilities. The outputs for the routine are: FX-the vector function (7x1), JAC-the Jacobian 
(7x5), and HES-the Hessian (5x5x7). Numerical results from the subroutine below have 
been compared with symbolically generated models for the same system using the Macsyma 
computer algebra system; no errors have been found. For this class of problems the user 
inputs data in the normal way, writes a math model for evaluating the function, and receives 
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the function, Jacobian, and Hessian without knowing any details about how the calculation is 
carried out. 
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VL2 Lagrange's Equations of Motion 

Lagrange's methods is a classical analytical dynamics method that requires the use of 
first and second order partial derivatives for the building the equations of motion for linked 
subsystems. A simplified analysis is presented because no rotating frame derivatives are 
considered. If m-generalized coordinates are available for analyzing the motions of the 
physical system, the velocity coordinates follow as: 

where q denotes the mxl vector of generalized coordinates and q denotes the time derivative 
of q. The velocity vector is used in building the kinetic energy, leading to: 



T = T(q 9 q) 

15 The calculations for Lagrange's equations are simplified by defining q and q as OCEA 

variables, and using these variables to evaluate the velocity and kinetic energy. The OCEA 
version of the kinetic energy becomes: 



20 



T, 


X 




~T 

,1,1 


T 1 

,1,1 








T ■ 
_ ,1,1 


T. . 

,1,1 _ 





(11) 



which is immediately useful for evaluating Lagrange's in the form: 

{T,,)-q = Q + T ,A T iM 

25 where the generalized force is given by Q q = ]T i*] v i q and the velocity partial derivative is 

obtained from the OCEA version of the velocity. This OCEA-based application produces an 
enormous reduction in the normal labor required for generating Lagrange's equations. 
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Similar labor reductions appear for Hamilton's method. A 2D example is presented in the 
next section. 



VL2.1 2D Example 

Assume that a cart of mass mi rolls without friction on a smooth surface. A linear 
spring restricts the cart motion in the x-axis direction. A massless rod of length 1 is attached 
to the cart and has a tip mass m 2 . The rod freely rotates about the attachment point to the 
cart. Gravity is the only assumed external force. The OCEA variables for the problem are 

(x,0, x,0) . The particle position vectors are given by: 



T o compute the kinetic energy the velocity variables are transformed to OCEA form, where 
the individual velocity components are given by 



1=(*,0) 

r 2 =(x + /*cos(0),-/*sin(0)) 



and the velocity variables are given by 

v,=(jc,0) 



v 2 = (x - 1 * sin(0) * 0, -I * cos(0) * 0) 
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V 2*-= 



x-l*sm{0)*0, 



-/*cos(0)*<9, 



0 

-l*cos(0)*0 
1 

-/*sin(0) 
0 

l*sm(0)*9 
0 

-/*cos(0) 



0 0 0 0 

0 /*sin(0)*0 0 -/*cos(0 

0 0 0 0 

0 -/*cos(0) 0 0 

0 0 0 0 

0 l*cos(0)*0 0 /*sin(0) 

0 0 0 0 

0 /*sin(0) 0 0 



The kinetic energy for this simple system is given by 

T = m t (v lx v„ + v iy v ly ) / 2 + m 2 (y 2x v 2x + v 2y v 2y ) / 2 
Using OCEA algebra leads to 



(/«, + m 2 ) x 2 1 2 - lm 2 sm(9)x0 + l 2 m 2 0 2 1 2, 
0 

-lm 2 cos(0)x0 
(m 2 +m i )x-lm 2 sin(0)0 
l 2 m 2 0-lm 2 sm(0)x 
0 0 0 0 

0 lm 2 sm(0)x0 -lm 2 cos(0)x0 -lm 2 cos(0);c 
0 -lm 2 cos(0)0 m 2 + m i -lm 2 sin(#) 
0 -lm 2 cos(0)x -/m 2 sin(0) l 2 m 2 
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where the vector and matrix partitions are easily obtained. In a real engineering application, 
the elements of T are numbers that are immediately used for evaluating the dynamic 
response. 
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VI.3 Equation of Motion Sensitivity. 

Given the first-order differential equation of the form 

M(x)x = f 

where the 3x3 matrix M(x) is given by 



lOOxj 5sin(jc,x 2 ) 50jc 2 jc 3 



M(x) = 

and the 3x1 vector f(x) is given by 



50xi 



SYM. 



X^ X 2 X^ 



^exp(-x 2 / 2) sin(jc 3 )^ 
-10cos(x 3 2 ) 



These equations are evaluated using OCEA algebra and the linear matrix equation is solved 
by Gaussian elimination. The results have been compared with the Macsyma computer 
algebra system and all first and second order partial derivatives have been correctly 
evaluated. Because of the limitations of space, only representative results are presented. 
Assume that (xi,X2,X3) = (1,2,3) as numerical values. The rate vector becomes 

'-0.34422£-02> 
jc= -0.28868£-01 
^ 0.17580£-02 , 

The partial derivative of the rate vector with respect to X2 follows as 

f 0.87449£-03 ^ 

0.28898£-01 
-0.189552^-02 



dx 

dX-y 



and the partial derivative of the rate vector with respect to X2 and X3 follows as 

f-0.36666^-02^ 



8 2 x 
dx 2 dx 3 



- 0.4773 \E- 02 
0.43682£-02 
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These results are typical. All of the first and second order partials have been validated by 
analytically inverting M(x), multiplying the result by f(x), and explicitly evaluating the 
partial derivative. 

VII Applications 

OCEA provides a rational process for generating sensitivity models by creating a new 
level of scientific and engineering software data abstraction and language extension. It has 
broad potential use for impacting the design and use of mathematical programming tools for 
applications in science and engineering. OCEA software replaces a time-consuming, error- 
prone, labor-intensive, and costly endeavor, with an automated tool that generates exact 
arbitrarily complex partial derivatives models. By developing the intrinsic operations at the 
scalar elemental level, the algorithm easily extends to Matrix applications. Future 
generalization will consider large-scale sparse applications. 
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